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ABSTRACT 


An existing code for calculating axial turbine performance 
using multiple stream surfaces was modified and made to run 
on the equivalent of an HP-1000 computer system. Calculations 
were made for the geometry of a 485 horsepower dual-discharge 
air-drive turbine for both on and off-design conditions. The 
results were compared with available data obtained at off- 
design speeds. Agreement of the flow rate and horsepower to 


within 5% was obtained. 
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I. INTRODUCTION 


A. DESCRIPTION OF THE TRANSONIC COMPRESSOR TEST RIG 

The Transonic Compressor Test Rig at the Turbopropulsion 
Laboratory (TPL) of the Naval Postgraduate School is shown 
schematically in Fig. 1 and consists of the following major 
components: 

1. Air drive turbine. 

2. Air supply system. 

3. Associated piping including throttling valves at the 
turbine and compressor inlets. 

4. Test compressor. 

ime drive turbine is a dual-flow axial air turbine with 
50% reaction. The geometry is given in Table 1. The profile 
shapes of the turbine rotor and of the stator blades are 
identical and the blades are of constant section along the 
radius as shown in Fig. 2. The stator has 31 blades while 
the rotor has 32 (to avoid resonant excitation from wake inter- 
ference). The two parallel stages of the turbine are designed 
for the following output and total inlet conditions: 

Pressure Ratio: 2.8 

Total Inlet Temperature: 640°R 

Flow rate: 10.85 LBM/SEC 


Horsepower: 485 HP 
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The compressor presently under test is a transonic single 
stage, axial flow compressor. It is instrumented for measure- 
ments of torque, mass flow rate, stagnation temperatures and 
pressures, case and hub wall pressures, and for unsteady 
pressure measurements in the flow field and at the walls. 

The Air Supply System incorporates an electric motor- 
driven multi-stage axial flow compressor manufactured by 
Allis-Chalmers. It can presently supply up to 12 lbs/sec of 
air at 3 atmospheres, at temperatures between 560°R and 660°R. 
The compressor is rated at 1250 HP and has a controlled 


variable speed drive. 


B. STATEMENT OF THE TASK 

The Transonic Compressor Test Rig was designed to pro- 
vide the means for obtaining experimental data in fundamental 
compressor phenomena. Following the present experiments, an 
experiment to investigate the onset of supersonic unstalled 
blade flutter is planned which would involve replacing at 
least the present compressor rotor by a rotating cascade of 
flat-plate blades. Such a rotor would not be able to produce 
the pressure ratios required to pump the required flow rates 
through the system. Therefore, it has been proposed, that 
a turbocharger compressor be fitted in series with the 
rotating cascade to provide the required flow through it. 
The turbocharger would also be driven using air from the 


Allis-Charmers air supply system. 
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In order to evaluate the feasability of the turbocharger 
installation, it is necessary to determine the mass flow rate 
required by the drive turbine to drive the test compressor 
at a given power and speed. The remaining air to drive the 
turbocharger turbine is then known and the selection of a 
commercially available turbocharger suitable for this 
application can be made. 

Thus, the performance of the air drive turbine must be 
known over the complete speed range. Of particular impor- 
tance, are the required mass flow rates for given values of 
horsepower. The problem, therefore, is to obtain the turbine 


performance map for all pressure ratios and speeds. 
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IIl. APPROACH 


A. BACKGROUND 

A search of the most recent literature revealed a number 
of analytical methods for the calculation of turbine off- 
design performance. The majority of these used in a finite 
element approach but little information on the relative suc- 
cess of these methods in practice was available. Two alternate 
methods, both used at the Turbopropulsion Laboratory, were 
those of M. H. Vavra and E. Macchi. Each was examined in 
detail. 

The method of Vavra, given in Ref. [1] is a one-dimen- 
Sional (meanline) approach using mathematical modelling and 
experimental data to express flow angles and losses. It is 
primarily a method to design turbine blading but may also be 
used to predict turbine performance for a given set of gas 
inlet and operating conditions when the blading geometries 
are specified. It is assumed that the axial velocity is 
constant along the blading from hub to tip. Vavra states 
that this assumtion is reasonable for blading in which the 
2p -LoO-nue. ratio 1S equal to or, less than 1.15. The ratio 
Potosi? and. 1.424 for the drive turbine stator and rotor 
blading respectively. It was thought therefore, that the 


method of Macchi might yield more accurate predictions. 
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Macchi's method is given in Ref. [2]. The method, imple- 
mented by Macchi in a computer program written for the IBM 360, 
was an extension of the work done by R. Eckert [Ref. 3] and 
R. Harrison [Ref. 4]. Eckert wrote a program, following a 
simplified three-dimensional analysis, which could be used to 
predict the performance of a single-stage axial flow turbine. 
Harrison improved the program by modifying the analysis to 
take into account streamline curvature. Both programs were 
based on the three-dimensional method developed by Vavra in 
Ref. [5]. Macchi's principle improvements to the program were 
to introduce the choice of various methods to calculate gas 
outlet angles and loss coefficients. Two methods of calcula- 
ting gas outlet angles are included; those of Ainley and 
Mathieson [Ref. 6] and Traupel [Ref. 7]. Five methods for 
eateculacine the loss coefficients can be selected; those due 
to Ainley and Mathieson [Ref. 6], Dunham and Came [Ref. 8], 
Baise [Ref. 9], Lonherr and Carter [Ref. 10] and Traupel 
met. 7). 

Macchi's computer program, as documented in Ref. ]2], was 
selected for performance predictions of the drive turbine. 

It should be noted that no card deck of the program was avail- 
able, and no results of using the program were available 


other than those included in Ref. [2]. 


B. ANALYSIS 

The method requires the following assumptions; 

1. There are an infinite number of blades in each blade 
row so that blades downstream do not affect upstream conditions. 


17 


= tmoiztime> 1 T2qu sostie. rom BB: mera Rito 2s 


"Yur + alam “ie 


— 
7 re 
-ofqmt:.,50rtSm sit ¢ 




















,0aF TAT ett 10? poretiw merger 
bre -[t ten) ttetod a, s00b 





“S10 out 26 Oe 
vitlvial PO? ,mETROTE 8 saomwe Suaitsd fh paris 

of boun sd bigos (hi  ehevignae Lene tensath- setds porninge i. 
-enidrey Wott Yates eee T gE 28 SDT MICT 19 T ad satel 
ae 2 BOR. SRI uUNEys Lbow yo mergota sdf bevorqmi noedrts 


amygorg dtoe .syimtevioo saiimessse tavoooe ornt Bee 


a TEV yt? Tee oh bofeem Lenetenemtii-eerdst ody ao boek 
aod 


Tow merged? 943 of atnomeyordm Siqtoniuq 2‘ insoaM 
Sloaleo of -yottewm epefyny To e2lons sit soubors i; 
-sivoiso to ebohsom owl .efsfenbitisoo e2e0f bas eeignas Jeive 
3 ho 30° 926: shyulont eye eslams sdefsuo 259 gee, 
' 9Ot ehodTSm-ave 5 . toh) fequsceT Bae [6 .202) coesiaae 
ere oft [astoeloe. 90 pac Seynadoratees seol sae snigeigogs 


_'& 365] eno bose mete? fd Stee) noeeer7aM bas vyelnia @ 


squat? fas [1 .368) 1eP'eo cam tteeegol gle - aan 


za ,..[ . $90 gf bornemsgb 26 ~maseerg radgmes 2*EdosaM 


-eo1d%ns Svievk edi Do Snolta rhage TREO TSS 765i. 5 ceroets 
-{feve sew wetyot0, sa 3a Bosh Bees Ga FRAT “beson od bTuode 
ida } [ow MEYSOTq Gat” earet ‘Yo 2s Tues T. or haw «St 


~ tS] 18% #1 Bobs ay Bi a2odt natty ss 
; r: 


earner ehataAn, A 


evOLtyamase smb wend fo ed? es sipet om iat . iT 

sbild Sone ab seboit PO. GE SEA, ne ene MT a 
aoe 
ak 


sm 


2. The flow is axisymmetric at locations where the 
equation of motion is solved. 

3. The flow is steady and adiabatic. Thus, the total 
enthalpy through the stator remains constant along a stream- 
line and the relative total enthalpy through the rotor remains 
constant along a streamline. 

4. All equations are solved at between blade row loca- 
tions. Increases in entropy occur in the blade row upstream 
of the stations where equations are solved and the entropy 
change along a streamline between blade rows is zero. 

5. The boundary layers on the turbine casing are not 
accounted for. 

The method of solution is as follows: 

1. Assume initial radial positions of the streamlines. 

2. Obtain the axial velocity distribution by solving 
the equation of motion at the stator outlet. The velocity 
distribution into the stator is assumed to be axial, and 
uniform 

3. Obtain stator loss coefficients. 

4. Check overall continuity and adjust the inlet Mach 
number aS necessary. 

S. Check the between-streamline continuity, and adjust 
streamline radial positions as necessary. 

6. Repeat this process for the rotor. 

7. Re-cycle all the above calculations, accounting for 


streamline curvature, and repeat until convergence is reached. 
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C. METHOD OF SOLUTION 

The computer code written by Macchi was originally run 
on the IBM 360 computer. The program consisted of a deck of 
over 2000 program cards plus over 60 data cards. Since the 
deck could not be located, it was necessary to re-type the 
program from the listing in Macchi's paper. However, since 
the IBM 360 computer was soon to be replaced in the period 
in which the work was to be carried out, an alternate 
computer was sought. 

The HP-1000 series mini-computer located at TPL was 
selected for two reasons. First, the machine used FORTRAN 
as did Macchi's program. Secondly, it would be a benefit 
to TPL to have the program immediately available on the 
laboratory computer. 

The first steps were to analyze Macchi's program, in 
detail, and then to run it using his example input/output. 
In analyzing the program it became obvious that the computer 
program listing given in Ref. 2, was not the one used to 
obtain the listed output. Numerous discrepancies were found 
in the listing, some of which would have prevented the pro- 
gram from running; others would have caused incorrect results 
to be obtained. A listing of these discrepancies is con- 
tained in Appendix E. When the program was understood and 
flow-charted, it was keyed-in at the HP-1000 computer ter- 


minal. However, modifications were required to accomodate 
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the program within the mini-computer disc-based operating 


system. 


D. MODIFICATION TO THE COMPUTER CODE 

Since there was no card reader, variable input data such 
as turbine speed had to be entered using data or specification 
Statements. This contributed in part to the most difficult 
problem, that of program size. The HP-1000 mini-computer uses 
a disc with a storage capability of 19.5 mega-bytes. However, 
the machine memory is only 124 K Bytes, of which only 29 K 
Bytes is available to a programmer. Also, the available 
memory is divided up, or partitioned into two 18 K and one 
11 K partitions, so that no single program can exceed 18 K. 
It was estimated that Macchi's program was over 100 K. So 
it was clear that the program would have to be modified if 
it were to run on the mini-computer. 

The first modification was to remove all subroutines 
from the program that were not actually used. It will be 
recalled that Macchi's program contained five methods for 
calculating loss coefficients and two methods for calculating 
gas outlet angles. It was decided that only the Traupel 
method of calculating loss coefficients would be retained. 
Traupel was selected for two reasons. Firstly, it was the 
method used by Macchi in his example calculations and there- 
fore the modifed program should still reproduce Macchi's 


results. Secondly, the method of Traupel is widely respected. 
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The method of calculating gas outlet angles was totally 
changed. Neither Ainley and Mathieson [Ref. 6] nor Traupel 
[Ref. 7] was used. Both methods required prohibitively 
large sections of computer code. The method selected was 
that of Vavra [Ref. 1]. 

Use of Vavra's method greatly simplified the program 
because this method predicts gas outlet angles independently 
of the inlet Mach number. Macchi's approach was to use 
Traupel's method which is dependent on the Mach number of 
the flow into the blade. 

The above simplifications reduced the program size from 
2257 lines to less than 1800 lines. However, this was still 
too large and the program could not be loaded without over- 
flowing the memory. 

The solution to the problem was found in program 
segmentation. In this process, the computer code is divided 
into a main program and several segments. Each segment is 
a "piece" of the original program. The segments are indi- 
vidually compiled and loaded. However, the segments are 
placed into memory only as they are needed to execute the 
overall program. Thus, a very large program can be made 
fo fun ifn the available 18 K partition. Since the present 
program was not originally intended for a mini-computer, 
segmentation was not straight forward. The method finally 


arrived at is detailed in Appendix C. Basically, the main 
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program consists of all the subroutines, while the three 
other segments contain coding which enables program flow to 
proceed in a logical manner. 

Successfully segmented, the program was run using Macchi's 
input. An output was obtained which agreed almost exactly 
with Macchi's results. All output quantities were within 1% 
of Macchi's quantities. The differences were, in all proba- 
bility, due to the different method of calculating gas outlet 
angles. 

After verifying Macchi's program, the drive turbine geo- 
metry was input and the program was run for a given set of 
operating conditions. The results are discussed in the follow- 
ing section. Note: The "verification" of Macchi's program 
amounted to verifying that the computer code now loaded into 
the HP-1000, was indeed Macchi's code. It was not known 
whether Macchi's output data were a good or bad prediction 


of performance since they were not compared with test results. 
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ITI. RESULTS OF AXIAL TURBINE PREDICTIONS 


A. USING BOTH COMPLETE AND MODIFIED PROGRAMS 

The drive turbine geometry was input and the following 
solution flow path was selected: 

1. Stator and rotor loss coefficients were functions of 
pressure ratio. 

ee the blockage factor, €*, used in the equation of con- 
tinuity was equal to the total loss coefficient. 

Four operating points were selected to test the validity 
of the program. Three were off-design points at which measured 
data were available and the fourth was. the design point itself. 
Table II contains details of the selected test points for 
Bun 

The program variables were then changed and the following 
new solution flow path was selected: 

ie stacor and rotor Loss coefficients were those calcu- 
lated by Traupel's method. 

7 the blockage factor, &*, was equal to the profile loss 
coefficient. 

After reviewing the results of Runs 1 and 2, a further 
modification was made to the program. The original program 
contained a subroutine which checked between-streamline con- 
fameity. if the total mass flow rate at the stator and rotor 


exits was not evenly divided between the five streamlines, 
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the radial positions of the streamlines were adjusted and all 
steps were recalculated using the new streamline positions. 
Hence, for Run 3, a subroutine was removed and the main pro- 
gram was modified so that between-streamline continuity was 


not examined. 


B. COMPARISON WITH MEASURED DATA 

tne nesults of Run 1, 2, and 3 are tabulated in Table T11. 
Run 1 showed predictions of mass flowrate which departed about 
6% from the measured data. However, the horsepower predictions 
were off by as much as 16.17%. Furthermore, the computer pro- 
gram was unable to reach a solution for the design point. 

Run 2 produced worse results as is evident from the table. 
Again, the program was unable to converge to a solution at 
the design point. 

Run 3 produced more acceptable data. Additionally, con- 
vergence to a solution was noticeably faster and a solution 
was obtained at the design point. Because of this, the method 
used in Run 3 was used to map the drive turbine performance. 
The computer program used to obtain the results of Run 3 is 
described in detail in Appendix A and is listed in Appendix G. 
The results of Run 3 are shown plotted in Figures 3 through 8. 

To obtain the plots in Figures 3 and 6, a value of the 
total inlet temperature was approximated by the method of 
VWavna as contained in Ref. [14].* It was assumed that the 


static turbine discharge temperature should not be less than 
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45°F (505°R). This corresponds to the approximate temperature 
at which condensation of moisture in the air, assuming 100% re- 


lative humidity, will occur. The inlet temperature was given by 


Total Inlet Temperature = Static Outlet Temperature 
la ea aes 

where n., the total-static turbine efficiency was assumed 
to be 81%, and Sp» Wes the total to static pressure ratio. The 
total inlet temperature corresponding to each pressure ratio 
os given in Table IV. 

The computer output corresponding to each point on Figures 
3 through 8 is contained in Appendix F. Only one side of the 
dual flow turbine was analyzed, thus, the resulting printed 
values of horsepower, referred horsepower, moment, referred 
moment, flow rate and referred flow rate must be doubled to 
obtain the actual turbine characteristics which have been 


plotted in Figures 3 through 8. 
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[V .. , DISCUSSION 


The agreement of both the predicted flow rateand the horse- 
power obtained in Run 3 with turbine test data was encouraging. 
It is to be noted however, that this agreement was obtained 
using a procedure which was conceptually incorrect. In Runs 
1 and 2, between-streamline continuity was checked and the 
streamlines were adjusted as necessary. In Run 3, between- 
streamline continuity was not checked, and as a result, the 
mass flow rate between streamlines was not precisely 25% of 
the total flowrate. It is noted however, that the deviations 
were less than 10.0% and while the radial positions of the 
streamlines varied by 10.%, the differences between predicted 
and measured output horsepower decreased from 24% to 4.5%. 
Since the enthalpy change on each streamline was computed 
using Euler's turbine equation, the total horsepower obtained 
by integration is sensitive to the streamline radial positions. 
On the other hand, the calculation of the overall mass flow: 
rate is primarily a function of the blade throat openings and 
inlet conditions of the flow. Consequently, in relaxing the 
requirement for between-streamline continuity, the output 
horsepower was changed significantly, while the overall 


flow rate was not. 
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Using this procedure, which preserves overall CONTINULEY , 
a performance map for the turbine was produced (Fig. 3-8) 
which agreed well with the off-design performance measurements 
made at lower speeds (Table III). It is noted however, that 
the inability of the program in its original form to predict 
the measured turbine performance is not explained, and both 
the program itself and the data input for the geometry should 
be closely re-examined. 

The difficulty in obtaining convergence to a solution at 
some operating points above the pressure ratio of 2.0 is 
likely to be the result of choking occurring on one or more 
of the streamlines. This was suspected but not fully explored. 

Finally, although the program was eventually made to run on 
the mini-computer, the time required to put the program into its 
final form was excessive since the original program was not 
written with segmentation in mind. When the segmented pro- 
gram was completed, only one operating point per run could be 
obtained. Thus, excessive time was spent compiling and loading 
the program. The execution time for the program averaged 2 
minutes at the lower pressure ratios and up to 30 minutes at 
the higher ones. This would be unacceptable if many points 


were to be examined. 
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V. CONCLUSIONS 


The program for calculating the performance of a single 
stage axial turbine reported by Macchi was revised, corrected 
and segmented and made to run on the Laboratory mini-computer. 
When applied to the geometry of the air-drive turbine of the 
Eempressor test rig, selecting specific options for the repre- 
sentation of loss coefficients, the revised program failed to 
converge when design-point test conditions where input. Also, 
the computed horsepower was in error by as much as 24% when 
the program predictions were compared with specific test data 
obtained from the rig at off-design (lower speed) conditions. 
The revised program did however closely reproduce the results 
given by Macchi in his original report for a specific turbine 
geometry. 

When the requirement that the computed stream surfaces 
be such that they divided the flow exactly into equal 25% 
increments was removed, the program converged satisfactorily 
for design point conditions and gave agreement with test data 
to within 5% in flow rate and horesepower at off-design 
conditions. 

The complete performance map for the air drive turbine 
was obtained with the program following this revision. Based 
on the favorable comparison with data so far obtained, the 
map is likely to describe the performance to better than a 10% 


uncertainty. This is considered to be satisfactory for 
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sizing the turbocharger for the proposed compressor rig 
modification. 

The following recommendations are made concerning further 

application or development of the computer program: 

1. The failure of the program to converge before the 
final revision was made should be analysed closely, 
and the final revision removed if possible. 

2. The geometrical input for the air drive turbine (which 
was taken from drawings) should be reexamined and the 
physical dimensions of the blade rows themselves should 
be measured. 

30) \Gonsideration should be given to putting the corrected 
original version of the program onto the IBM 370 com- 
puter so that, when successfully operating, a turbine 


map can be calculated with a single load. 
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TABLES I 


TURBINE GEOMETRY 


(see Figure 2; Dimensions in inches) 





























STATOR 

Hub Radius 20164 
Mean Radius 5.98: 
Tip Radius S027 - 
Blade Chord 1.003 
Blade Suction Side Radius of Curvature 2.8065 
Maximum Blade Thickness SLZSZ 
7.&. Projected Thickness 03 
T.E. Normal Thickness .0186 
ROTOR 

Hub Radius 2609'S 
Mean Radius 5.265 
Tip Radius 5.857 
Blade Chord 1.003 
Blade’Suction Side Radius of Curvature 2.8065 
Maximum Blade Thickness ee 
Z-—. Projected Thickness 03 
Toe. Normal Thickness .0186 
Tip Clearance .Ol(estimated) 
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TABLE TITEL 
MEASURED/DESIGN DATA USED TO VERIFY THE PROGRAM 


LBM 


POINT BPM Tin p)lour (ry to cpst)p.r. Sec) H-P- 


1 Peord” “579°. Z S50.6 25.90 1.602 soo4e 210.1 
2 L5200~.-.557..-4 917-8: — 20.43 1-.390 42-698 —63-27 
3 £13004 578.9 306.8 027.13 1.846 72033 19230 
4x 30500 640.0 ae 41.16 Lav8 10.85 485 


*Design Point 
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TABLE III 


COMPARISON OF PREDICTED TURBINE PERFORMANCE 
VS_ MEASURED PERFORMANCE 























RUN I 
FLOWRATE HORSEPOWER 
POINT PREDICT. MEAS. ‘DIFF. PREDICT: MEAS *DIFF. 
1 3. o0 5. DAZ 6.09 2365 ETO 9.05 
2 4.74 4.698 0.89 Sag) 63e27 bcd Caf 
3 7.04 7.0355 0.009 163.64 72.0 4.86 
4 N.C. L055 25 NaC: 485 === 
RUN 2 
i 6.06 5. 542 S205 90.92 EO 17.4 
2 4.90 4.698 4.29 49.76 OSiac? @ aia S5 
3 7.30 PsOSS 3.80 1350576 2) Let 
- Nw 10..85 —— N.C 485 aaa 
RUN 3 
1 5, 82 5.542 5.02 DS 752 LINO Gat 2.74 
2 4.66 4.698 0.81 61. 916 65527 2109 
B 7.04 72055 OL LO 7-9 S68 LT ZO 4.47 
1 10.40 10.85 4.15 444.18 485 8.42 


NC: Computer program would not converge to a solution after 
a large number of iterations. 
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TABLE IV 


VALUES OF ASSUMED TOTAL INLET TEMPERATURE FOR EACH 
PRESSURE -RATTOSGIVEN IN FIGS. 3, 5, 6,..AND27 


PRESSURE RATIO TOTAL INLET TEMPERATURE (°R) 
1.4 545.5 
i.0 502.0 
1.8 Sil 2S 
Z.0 592.0 
D2 603.6 
Zt OA 5 33 
Zo GIZiO. L 
Le 636.6 
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FIGURE 1: SCHEMATIC QF THE COMPRESSOR TEST RIG, WITH 
PROPOSED MODIFICATIONS 
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APPENDIX: A 


DESCRIPTION OF THE COMPUTER PROGRAM 


A-1. INTRODUCTION 

To enable the program to run on the laboratory computer, 
the program was divided into 4 parts; a main program and 3 
segments. A detailed discussion of program segmentation on 
the HP-1000 computer series is contained in Appendix C. In 
the description which follows, the program is treated as if 
it were one large program with many subroutines. 

The description follows the individual steps from 
start to finish in the analysis. A program flowchart is 
given in Figure A-1 and the FORTRAN symbols used in the 
program are listed in Tables A-I to A-IX. 
A-2. - DESCRIPTION 

A=2 A> Input. Data 

There are 4 basic categories of input data; tur- 

bine geometry, operating conditions, special data and pro- 
gram control parameters. Since there was no card reader 
input device on the computer, all datawere entered using 
either data or specification statements. Explanations of 
the turbine geometry, operating conditions, special data 
and program control parameters are found in Table A-I through 
A-V. The nomenclature for the blading is given in Figure 
A-2. 
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A-2.2 Initial Geometric Calculations 

The first calculation performed is to establish 
the 5 streamline locations at the stator inlet (station 0). 
The streamlines are initially positioned such that there are 
equal areas (25% of the total flow area) between them. Next, 
blade heights of the stator and rotor are calculated using 
the hub and tip radii of each blade. Blade spacings for the 
stator and rotor are computed at 3 streamlines; hub, mean 
and tip. The blade spacing on the mean streamline for the 
stator is given by 
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S = Z. Rm (A-1) 
where S = Blade spacing 
Z = Number of stator blades 


Rm = Mean stator radius 
Ao2. 56> Caleulation of “Gas Outlet Angles 
Subroutine VAVRA calculates gas outlet angles 
for both stator and rotor. The method is that of M.H. Vavra 
[Ref. 1]. The equation programmed in the subroutine is 


Oe bared uy a (|- 2) oe 





where a = Gas outlet angle 
a = Throat opening 
S = Blade spacing 
t,= Projected trailing edge thickness 


This method is much simpler than that used by Macchi since 
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there is no variation in outlet angle with Mach number (for 
sub-sonic conditions). Therefore, once calculated, the stator 
and rotor exit angles remain unchanged. Subroutine VAVRA 
computes exit angles for the hub, mean and tip streamlines. 
The outlet angles at streamlines two and four are computed 
later in the subroutines STATR and ROTO2. 

Before printing the input data, the program 
calculates the mean throat opening for the stator and for the 
rotor. The ten equally spaced radii and corresponding throat 
openings (part of the input geometry) are fitted with a fourth 
order Chebyschev polynomial. A throat opening corresponding 
to the mean radius is thus obtained. In the present appli- 
cation of the program to the drive turbine, the mean throat 
opening was obtained from the design drawing of the blading 
shown in Figure 2. It was assumed that the throat opening 
varied linearly with radial position and hence the throat 
openings at other radii could be calculated. The resulting 
throat openings are shown in the computer output under the 
neading sof “Input Prints". The design values, of the stator 
and rotor throat areas were obtained from the original 
design notes of M.H. Vavra. 

A=2.4% Galculation of theFlowRate 

Subroutine CHAN is called to calculate the mass 

flow rate entering the stator. The equations used are as 


follows: 
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TS THe 
} + Y-l Mo” 
Mise 4 Oe: ty Pods (A-4) 
‘i 
ae are ear (A-5) 
Ute a ees 
pate Ry ae 
2 2 
A= 17 LR wwe — RvB C7) 
Dn = Gee Pal, (A-8) 
2 | Re (A-9) 


YY 
E = 
BER Terps 
is the reference (dimensionless) flowrate 


m 
nee 
and is used to check overall continuity later in the program. 
A-2.5 Solution of the Equation of Motion for the Stator 
Subroutine STATR is called to solve the equation 
The equation of 


of motion for the stator outlet conditions. 


motion which is programmed is as follows: 
Beye) costa -(H AT m 32). (tee 
C, Heosol A shivrol ee 
AY, (A-10) 


C, Cosel aH ‘z 
ae Van AX ‘i Var 


45 









- i 1, - 
; d 
f P. 
{ gi } 
# — 
& 
n\ g > rs vd 
re l 
- a wi ao 
— Ty eA 
- ) 
| AS = Mt 
— -~ a Rance Sd 2 e 
a ry 
my ” ire = 2 vt i 
Se ee eee ‘ - a el 
ae } t + " 2 - 
, e . 
- (erel 7 1 bi) . WIA : - ry s ‘3 a * 1 . 
. —_ F eke = a ; as ter” 
a g od-m! 20051 vetuntsmos ifavrevo Anata 62 bee ae: 


27870 Sas 304 noisem 1o fo lpaas ose 19, ho ttut oe 2,8 
ofiiops 383 eylose ct bheliaos =) BATE wLEsuOTdUR a 


te Abvtavgpe sa bets REGS seliuc Tosede edt rot 20k 


real iat hid ea nara dokdw 





where C, = A Je J (a constant to convert H, the enthalpy 
from BTU es 


LBM TO sec 


VV ) = Axial velocity at a streamline 
a (3 Axial velocity at mean streamline 





x = RC) . Streamline radius 
Ee ine Mean streamline radius 
Won Vee 
ne et = st coset 
= — dn Y, 7 V, z 
aX, cl X = ee 


Lee , eer 
C, H costa (1- & ) 
6-= Stator loss coctficient 
(which is initially assigned an 
estimated value) 
The derivation of this form of the equation of motion is 
given in Appendix B. However, at this stage of the analysis, 


the streamline curvature is assumed to be zero. Therefore, 


the equation of motion becomes: 


2S ae ¥ + 7 peras oleh : wae, a C, cos‘, 


SIN of; + 
dX AX =X Yoo 
Bee by Cr Hicos os ds* | 
me |! Wee or se (A-11) 


The equation of motion is solved when the value 
of Y; at each streamline satisfies the equation. The solu- 
tion is to first put the equation in the form: 

a 
Y, 
d (Ln%") alee 
dX, 


where I(X) consists of the right hand side of equation (A-11). 


(A-12) 


Integrating equation (A-12) yields; 
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Xo 
where &£nc* is the constant of integration when x = 1 and Y.= 
1 


1. With these boundary conditions Eq. (A-13) gives 
1 


Ln c* =. > ( t(xjid Xi (A-14) 
XO 
msing Eq. (A-14) in Eq. (A-13), 


eye 


Xo 


( | 
IOdX,-| LOdd X, (A-15) 


which becomes . 
t 


Ln Yee ready we 


Taking the inverse natural log and the square root of both 


= xX 
sides 4 ( 20odn 
a | 
ae (A-17) 


Equation (A-17) is the form of the equation of motion solved 
in subroutine STATR. Solution of the equation gives five 
values of Ls and thus the value of the axial velocity at 
each of the five streamlines. Initially, the value of Y, 
is taken to be 1 and the value of — is taken to be zero. 
In succeeding iterations, the eallemineed value of ae is used 
to obtain a new value of — and so on. 

After Sofie: oda five values of te: the ‘stator 


exit conditions are calculated at each streamline from the 


geometry of the velocity diagram. The convention for positive 
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and negative angles and velocities is defined in Figure A-3. 


The required relations are the following: 


[Un GUS is (A-18) 
Vin, = Va, > TAN, (4-19) 
Vi = Vay J 60s &: (A-20) 
Vee -Va, | AR /2L | (A-21) 


Where L is the axial distance between stations and AR is the 
change in radial position of the streamline. 2 the radial 
component of velocity, is taken to be zero at this stage in 


the calculation. 


Vi = vis Ve 
2. (A-22) 
eee Vi 
Tus = Tro - roe Ts 
L (A-24) 
| §s 
A, a if Pro Y (8225) 
os P Tits r-{ 
Pi = TO ! a (A-26) 


M, = Vi /4¥qe RT ae 


After the above quantities have been calculated at each 


streamline, subroutine STATR returns to the main program. 
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A-2.6 Calculation of the Stator Loss Coefficients 
The calculation of the stator loss coefficients 
at each streamline is accomplished by subroutine ALOS1. 
The method of solution to obtain these loss co- 
efficients is that formulated by Traupel [Ref. 7]. In 
Traupel's method, the value of the total loss coefficient is 


given by 
E ae +E ame (A-28) 
total profile wall remaining 


une calculation of € requires 9 subroutines. Figure A-4 


total 
describes the connection between the subroutines and subroutine 


ALOS1. 

The first step is to obtain the value of the 
total profile loss coefficient, EDt ED is defined by Traupel 
to be 


5, ; ee ms ais z t, (A-29) 


where Epo = initial value of the profile loss coefficient 


a mach number correction factor 
ais trailing edge thickness correction factor 
a = loss coefficient due to mixing losses and separa- 


tion losses 
= loss coefficient due to fan losses 

The total profile loss coefficient is calculated 
am the following manner. First, data for initial profile loss 
(E50) as a function of gas outlet angle (a,) for various values 


of gas inlet angle (a,) is read from an array (Fig. A-5). 
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This is done by subroutine TRAU1 and functions XPO and YC. 
The values of EDO are contained in two arrays XPOl (5, 8) and 
XPO2 (6, 8). This is because the data shown plotted in Fig. 
A-S has been divided into two sets. One set is for values 

of a, between 40° and 80°. The other is for values of a, 
between 80° and 170°. The FORTRAN symbols for the two ranges 
of values of a; are ALFOI1(I) and ALFO2(I) respectively. The 
FORTRAN symbol for the gas inlet angle is ALFl (J) once the 
data pints selected from the plot are entered, fifth and sixth 
degree polynomials respectively are fitted through the data 
points. The value of Eno can then be determined for given 
values of a, and ay. 

The mach number correction, Xn is obtained from 
Fig. A-5. Subroutine CSIM calculates the value of Xn using 
straight line approximations of the plot. 

Subroutine CID calculates the remaining terms in 
the expression for EDe These are Xo; ae ee: They are ob- 
tained from the data in Fig. A-6 using the linear interpolation. 
The abscissa of the curves for X. and En is either’ f prii-£ 


where f is defined as S 


f abs Sick, (A-30) 


= 
Io 
oO 
8 
@ 
O 
i) 


normal trailing edge thickness. 


ct 
i} 


blade spacing. 
@y = gas, outlet. angle. 
The Joes, coctticient due to wall friction, € , 


W 


is calculated using 
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é = ¢ ™e T SINK 

Ww Po P ZL (A-31) 

where t - blade maximum thickness 
£ = blade height 
This equation is programmed in subroutine CSIW. 
The value of ER is obtained using subroutine CSIR. 

ER is defined by Traupel to be an all-inclusive loss coeffi- 
cient which accounts for any remaining losses not previously 


defined. It is written as 


= X, tee (A-32) 


ERO is an initial value of ER which depends on the value of 


®?, where $¢ is given by 
Vi SIN, 
p = VU (A-33) 


in which v, = true velocity of gas 


Vv blade speed 
A plot of ERO vs ¢ is shown in Fig. A-7. The correction X 


is a function of s/2 where 


Ss chord length 


R 


blade height 

aiie is Obtained using the data inythe lower half of Fig. A-7. 
The total stator loss coefficient is computed 

for 3 streamlines; those at the hub, mean:. and tip. 

The loss coefficients at streamlines Z and 4 are 


obtained by linear interpolation. 
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A refinement to the stator loss coefficient may 
be applied depending on the input value of one program control 


parameter. The following 3 variations of bs are available: 
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where € = loss coefficient calculated using the method of 
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The values of the program control paramenter 


required to select between options are given in Table A-V. 
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Before returning to the main program, subroutine 
Poeat Gabeulates a value of €* which is a blockage factor to 
be used in the equation of continuity. There are three ways 


to define &*; they are as follows: 


hek.s 
So (A- 38) 


(A359) 


f = be (A-40) 


A-2.7 Solution of the Continuity Equation After Re- 
turning to the Main Program 
The overall continuity at the stator exit. 1s 
checked. Subroutine FLOWR performs this task. The flow chart 
for FLOWR is given in Fig. A-8. In FLOWR the mass llawerete 
required by continuity is checked against the calculated mass 
flow rate. If the calculated flow ratedoes not agree with 
that required by continuity, adjustments are made to the axial 
velocity and/or the inlet Mach number, as will be explained. 


The mass flow required by continuity is 





“package f M REF 
Meeap = SAAR 
ee ate (A-41) 
where i = reference mass flowrate as computed in subroutine 


CHAN 
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Ss 
a = mean stator throat opening 
aa = mean stator radius 


The mass flow rate at each streamline com- 


puted in this subroutine is 
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ACT 
where Z is an area reduction coefficient defined by 
H e KK 
ral +x re 
ae f 
H | (A-43) 


Z gives the percentage of flow area between the blades over- 
which it is permissable to assume a uniform velocity. The 
boundary layer on both sides of the flow limits the available 
flow area and the backage facotr, Z. accounts for this. 
Equation A-43, Z is seen to be afunction of the energy para- 
meter H*** and —€*. &* is the value of the loss coefficient 
returned from subroutine ALOS1. The energy parameter is 


defined as 
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63 Ener thickness 
keEF = cia = = 
: Displacement thickness (A-44) 
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H*** can be written as 
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(A-45) 
where: 
m= .15 
Xo = ee a for unchoked flow 
E Pan 
To 
y-1 
x = TP iT! x for choked flow 
and : Y 
= Wik Star: 
Perit ~ | yi] 


The derivation of Z and H*** os given in Appendix B. 
The expression for ¢, the flow function, 


for unchoked flow is 


(A-46) 





and for choked flow is 
ot ica 
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After calculating for each streamline, the flow 
rate is integrated from hub to tip and the resulting value is 
compared with ere If the two values of flow rate agree to 
within a specified tolerance (see Table A-IV) continuity is 
considered to be satisfied. Then, after calculating the total 
percentage of mass flow between adjacent streamlines, sub- 
routine FLOWR returns to the main program. 

If the flow rates are not within tolerance the 
program checks to see if the actual mass flow is to high. If 
it is to high, the value of the axial velocity is lowered 
proportionally to the difference between the actual and 
required flow rates. 

If the actual flow rate is too low, the procedure 
is more complicated. First, the flow is checked to determine 
whether choking has occurred. Streamlines one and five are 
checked. If the flow is in fact choked at those streamlines, 
the inlet Mach number is lowered and the program loops back 
to recompute the reference mass flow rate and repeat the 
complete procedure. 

If the flow,is not choked, the axial velocity is 
raised proportionally to the difference between actual and 
required flow rates and subroutine FLOWR returns to the main 
program. 

A-2.6 Calculation of the Rotor—Iniet Conditions 

Continuity having been satisfied through the 
stator, the rotor relative inlet conditions are calculated. 


In subroutine ROTO1, the following expressions are used: 
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He = heya 


Where Trp? PrE and He are equivalent temperature, pressure 


(A-56) 


and enthalpy respectively. 
A-2.9 Calculation of the Rotor Exit Conditions 
Calculation of the rotor exit properties follows 
the same procedure as was used to compute the stator exit 
properties. The process is outlined here with notable 
aisrerences explained. Subroutime ROTO2Z calculates the 
Eecor Exit properties. A flowchart of ROTOZ is’ given in 
Fis. A-9. 
The first step in ROTO2 is to solve the equation 
of motion for each streamline. The equation of motion in 


terms of relative quantities is 


Ya" Vas 


Bie ens Bs. " Gs cos Bz dHe 
Ya Va Van . ate 


i 6: He cos’ Ba : ee (A-57) 
Ya? Vas d Xz 
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At this point in the calculation, streamline curvature is 


neglected. Hence, Eq. (A-57) reduces to 


a Ys ) ios ol 62 X 
A Xa : sede 


a = 
dX2 Xo ae Ba 





4 Um Cos Ba SIN a iz ey Sea U3 Cos Ba 
Ya Vey ea : 


C.cos"B2  olHe Ae C. He cos*b2| Se 
ee, aX2 iG ae 1X 


(A-58) 


The derivation of Eq. (A-57) is contained in Appendix B. 
Equation (A-55) is similar in form to Eq. (A-10). Hence, the 
method of solution is identical to that employed by subroutine 
STATR. However, after solving the equation, the value of oe 
at each streamline is examined to determine whether or not LE 
falls into the range ou ae Values of = greater than 
2.0 are set equal to 2.0 while those less than .2 are Set 
equal to .2. Successive values of Lis at each streamline are 
compared, and when the values of successive iterations are 
within a specified tolerance (see Table A-IV), the iteration 
ends. The values of =e are used to calculate the rotor exit 


conditions using the following equations: 


Vax = Vaa (3) 1é (A-59) 
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(A-60) 


(A-61) 


(A-62) 


(A-63) 


(A-64) 


(A-65) 


(A-66) 


Subroutine ROTOZ, then returns to the main pro- 


gram. 


After calculating the rotor outlet conditions, 


the rotor loss coefficients are computed. 
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Subroutine ALOS2 
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calculates the rotor loss coefficients following the process 
used in subroutine ALOS1 for the stator losses. The principle 
exception is that a tip clearance. loss.is also.calculated and 
added to the total loss coefficient. The tip clearance loss 
coefficient is obtained from subroutine ALEAK which uses a 
straight line approximation to the curve shown in Fig. A-10. 
Subroutine ALOS2 also computes values of &* and one of the 
three refinements to ER- 

Subroutine FLOWR is called to check continuity 
mencne rotor exit.. If continuity is_satisfied, the program 
continues. If not, the same procedure is followed as pre- 
viously explained for the stator outlet (Fig. A-1). 

A-2.10 Accounting for Streamline Curvature 

All calculations to this point have neglected 
streamline curvature and assumed that the streamlines remain 
fixed through the stator and rotor (Fig. A-11). The radial 
shift in a streamline between stator inlet and rotor outlet 


can be written as 


ra Rca Ts RoToOR (A-67) 


INLET OUTLET 
This is the net radial shift in a streamline between stations 
a (sfator inlet) and "2" (rotor outlet). It is shown in 
Section 16.4 of Ref. [5] that the radial shift in a streamline 
between the stator and the rotor (station 1) can be written 


as 
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OvuTLET INLET OUTLE1 


The angle between the meridional velocity Ue and the axial 


velocity y is dX. The radial velocity es can be expressed as 








(A-69) 
“eee erom Fig. 16(1) of Ref. [5], it follows that 
. = AR 
TAN "a al (A-70) 
Mermenckq. (A-68) in Eq. (A-67), 
V Ele 
Ve = 7 Ve aie (7) 
where on = Average streamline slope 
Also, from using Eq. (A-68) 
gees 
= (A-72 
Cos ) 5 - 
| AR™ +(@L) 
Rearranging; 
(A-73) 
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The remaining term used in the calculation of streamline cur- 


vature (Section 16-4 of Ref. [5]) is 


6R 
ie 


where K is the so called curvature factor. It usually has 

a value between 4 and 6 and in the program its value is taken 
to be 5. Having calculated cosd, AR and 6R, the program 
repeats the solution process. However, the only quantity 


which is unchanged is the reference mass flow rate m In 


Tels 
subroutine STATR the equation of motion is solved, this time 
accounting for streamline curvature. The same is true in 
subroutine ROTO2. 

The flow path of the program is identical to the 
section which did not account for streamline curvature. Next, 


the program computes an average pressure ratio at the rotor 


outlet using the expression 











Po, = P “ on P, 
Pro Fro | STREAMLINE 4 Pro SL. 
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P, Pa, Pa 
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(A-74) 


He. this pressure ratio is withinga, specified tolerance to 


the actual pressure ratio (which is input data) the program 
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proceeds to the final stage of the calculations. If not, the 
inlet mach number is adjusted by an amount which depends on 
the difference between the calculated and specified peesane 
Bataos. if the calculated pressure ratio is too high, the 


Mach number is lowered using 


te » Bressune Ratio Ditference 
Mo MS 18 (A-75) 


If the computed pressure ratio is too low, the Mach number 
is raised using 


M =m + Pressure Ratio Difference (A-76) 


Oo O 


In both cases, the program loops back to subroutine CHAN and 
proceeds to compute a new reference mass flow rate based on 
the new value of the inlet Mach number. The entire process 
is then repeated until the pressure ratios agree within the 
specified tolerance. 
A-Z2.1d, Final Calculations 
Staeor ana rotoryoutlet] conditions not previously 


calculated are computed as follows: 
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(A- 86) 


(A- 87) 


(A-88) 


(A- 89) 
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The turbine horsepower is obtained by integration. 
The Ah term at each streamline is weighted by the percentage 
of mass flow at that streamline. The product is then inte- 
grated from hub to tip and result, An, as used in the tur- 


bine horsepower equation 
Ah] - Mm 


550 (A-95) 


eee. = 
The moment is calculated using 


(H.P, )( 550) 


OW) (A-96) 


Referred horsepower, moment, mass flowand RPM are calculated 


using 
H.P. 


H.F, REF = Oo S (A-97) 


M ea (A-98) 
™ ger = SS (A-99) 


ReMeee- > (A-100) 
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The values of the total-static efficiency, total- 
fatal Efficiency, total-static pressure ratio, total-total 
pressure ratio, head coefficient, blade/jet ratio, r* and 
inlet mach number are then averaged. 

With ail calculations completed, the results are 
printed under the heading "STATOR SOLUTION", "ROTOR SOLUTION", 
and "OVERALL TURBINE CHARACTERISTICS". 
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TABLE |. A-I 


TURBINE GEOMETRIC INPUT DATA (STATOR) 
(see Figure A-2; Dimensions in inches) 


FORTRAN SYMBOL DESCRIPTION 

ZS Number of blades 

RS (1) Hub radius at stator outlet 

RS (3) Mean radius at stator outlet 
RS (5) Tip radius at stator outlet 

2 Blade chord (mean) 

CI Blade chord (hub) 

co Blade chord (tip) 

E Blade curvature (mean) 

EI Blade curvature (hub) 

EO Blade curvature (tip) 

T Maximum blade thickness (mean) 
TI Maximum blade thickness (hub) 
TO Maximum blade thickness (tip) 
TE Projected T.E. thickness (mean) 
TEI Projected T.E. thickness (hub) 
TEO Projected T.E. thickness (tip) 
TN Normal T.E. thickness (mean) 
TNI Normal T.E. thickness (hub) 
TNO Normal T.E. thickness (tip) 
Al(1-10) Ten values of throat diameter at 


10 equally spaced Yadia 
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FORTRAN SYMBOL DESCRIPTION 


AL Blade camber line length (mean) 
ALI Blade camber line length (hub) 
ALO Blade camber line length (tip) 
BC(1) Hub radius at stator inlet 

RC (3) Mean radius at stator inlet 
RC(5) Tip -Tradius~at-stator-iniet 
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TABLE *A=Ti 


TURBINE GEOMETRIC INPUT DATA (ROTOR) 
(see Figure A-2; Dimensions in inches) 


FORTRAN SYMBOL DESCRIPTION 


ZR Number of blades 

RR(1) Hub radius 

RR(3) Mean radius 

RR(5) Tip radius 

CR Blade chord (mean) 

CIR Blade chord (hub) 

COR Blade chord (tip) 

ER Blade curvature (mean) 

EIR Blade curvature (hub) 

EOR Blade curvature (tip) 

TR Maximum blade thickness (mean) 
TIR Maximum blade thickness (hub) 
TOR Maximum blade thickness (tip) 
TER Projected T.E. thickness (mean) 
TEIR Projected T.E. thickness (hub) 
TEOR Projected T.E. thickness (tip) 
TNR Normal T.E. thickness (mean) 
TNIR Normal T.E. thickness (hub) 
TNOR Normal T.E. thickness (tip) 
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FORTRAN SYMBOL 


A2(1-10) 


ALR 
ALIR 
ALOR 
CV 
CK 


DESCRIPTION 


10 values of throat diameter at 
10 equally spaced radii 


Blade camber line length (mean) 
Blade camber line length (hub) 
Blade camber line length (tip) 
Axial distance between stations 


Curvature Factor 
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TABLE <A-III 


TURBINE OPERATING CONDITIONS (INPUT DATA) 


FORTRAN SYMBOL DESCRIPTION 
AMC Assumed inlet Mach number 
Assumed stator exeit Mach number 
AMS (absolute) 
Assumed stator exit Mach number 
AMR (relative) 
PTO Total inlet pressure (PT0) 
TTO Total inlet temperature (To) 
PR Total-static pressure ratio 
RPM Operating speed (RPM) 
VA1(3) Assumed axial velocity in stator 
VA2(3 Assumed axial velocity in rotor 
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TABLE " A-IV 


SPECIAL INPUT DATA 


FORTRAN SYMBOL DESCRIPTION 
Toler mce for convergence of 
TOL 1 : Leah: 
equation of continuity 
TOL 2 Tolerance for between-S.L. 
continuity (not used) 
Tolerance in pressure ratio 
oes convergence 
TOL 4 Tolerance in equation of motion 
convergence 
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TABLE A-V 


PROGRAM CONTROL PARAMETERS 





FORTRAN. SYMBOL POSSIBLE VALUE EFFECT /MEANING 
IND 1 routines STATR, FLOUR, 
ROTO2 
1 No printing in the above 
ICL 1 Rotor is shrouded 
1 Rotor not shrouded 
ICOZ 1 E = Eo 
6 & = € (Y Pressure Ratio) 
ae: 
ICON 1 E= SS ror ar 
Z = EpROFILE 
$ © = SToTAL 
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TABLE | A-VI 


FORTRAN SYMBOLS IN THE MAIN PROGRAM 

















FORTRAN SYMBOLS DESCRIPTION 
BESP gee pient -(.8) 274 
OI Stator throat opening (hub) 
00 Stator throat opening (tip) 
OIR Rotor throat opening (hub) 
OOR Rotor throat opening (tip) 
0 Stator throat opening (mean) 
OR Rotor throat opening (mean) 
ANG2I Stator gas outlet angle (hub) 
ANG20 Stator gas outlet angle (tip) 
BETAI Rotor gas outlet angle (hub) 
BETAZ Rotor gas outlet angle (tip) 
G Grav. constant, o2 7174 Peel 
GJ 778-16" FT sLBF/ BTU 
EXP1 se i eee 
EXP2 y-Ty 
ERRE Gas constant, 53.3459 oer 
EMME Molecular mass, 28.970 LBM/LB MOLE 
GAM y, Ratio of specific heats 
ETAT Total-total efficiency 

ETAI Total-static efficiency 

ETAS Stator blade efficiency 
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FORTRAN SYMBOL DESCRIPTION 





























ETAR Rotor blade efficiency 

RSTAR Theoretical degree of reaction 
ALOS Head coefficient 

BLAJE Blade/jet ratio 

DR1 Radial shift of steamlines 

AMW1 Stator exit relative Mach Number 
AMS1 Stator exit absolute Mach Number 
AMV 2 Rotor exit absolute Mach Number 
AMR2 Rotor exit relative Mach Number 
DELH A 

HP Horsepower 

AMOM Moment 

THETA 8 

DELTA 6 

HP1 Referred H.P. 

AMOM1 Referred moment 

RPM1 Referred RPM 

WLBM1 Referred mass flow rate 

ETAS Average total-static efficiency 
BETA6 Average total-total pressure ratio 
ETA6 Average total-total efficiency 
AKISS Average head coefficient 

RSTARS Average theoretical degree of reaction 
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TABLE A-VII 


FORTRAN SYMBOLS IN SUBROUTINE CHAN 





FORTRAN SYMBOLS DESCRIPTION 

TTO Tro? total temp. at station @ 
AMC Inlet Mach number 

PTO Pro» total pressure at station § 
RC (I) Streamline radii at station § 
WLBM M, required mass flow, pAV 

pa Static temperature 

Ve Velocity 

PG Static pressure 

RHO Gi5.censity of air 

WCHAN sae reference mass flow 
WPERO % of M at each streamline 
































f 
YAO AT Ti wm > 
AAD BY ITOK ME Bd 
; ; VOI (ea290 
- - ome _ ibe ea ll sa ~ en ‘ 
— mE NT a TS RNR RUN Ls , i A “ 
cofiet: 0 .quee fntns Feypt 
Se ee “ 
Tsdiwk aah sadset 
—_ pane lge ve - ith neem a cat geoe —— : _ 
/ tolsete Th Sstweetsyq feso? 
. ... . i i * oT 
~~ my ES Ef 26 te a A AR go = y me 
) taktete gs Feber sal lagorke 
a —— os = = ee So — 
i, F ' a ae : — ¥ 
-- Vat WoL? ebSagi HetLopsay iM 
as ams lia et ee | mmm a a ee Ra 
8 | everett at page 
— ene, Soe diee——syyepenssemmdsanssibp i lchae a rage 
ytioolsyv 
— — — Se — NT CR 
Sivezgoets 2itewe 
a a A . - : Se mera 
118 40 yri2peh a 
—. —— eee errhp ——— ne ee ; 
i, , _ - ° ve ‘ 
"OL? 2EGh Sone sete9 gage 
1 Pa aa pm aie i ae 
og ate C4 ao a aa 
SNC LMESITIS2 PFoas ge M 345 @ 
——— A | TTT TERRI _— eer . 





TABLE A-VIII 


FORTRAN SYMBOLS IN SUBROUTINE STATR 











FORTRAN SYMBOL DE SCRI PTION 

ALFA1 Stator gas outlet angle 

x Ratio of streamline radius/ mean 
radius 

AMS Mach Number at station l 

- Static temperature 

P Static pressure 

V1 Absolute velocity 

VAl Axial velocity 


Ratio of axial velocity to mean 


axial velocity 
S Entropy 


DSDX Entropy gradient between streamlines 


vul Tangential velocity 
PRAT (Total-static pressure ratio) 
£i1S Tits 


DALF 7a 








RSF Mean stator radius 

R : R 
DELR Stator .in = rotor out 
ZETAPS § 

Pp 
ZETAS § 

S 
VR1 Radial velocity 
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TABLE A-IX 


FORTRAN SYMBOLS IN SUBROUTINE TRAU2 

















FORTRAN SYMBOL DESCRIPTION 
CSIP Xp> correction to 
R Se initial profile loss coefficient 
ANG1 Gas outlet angle 
ANG2 Gas inlet angle 
Rl Xa? Mach No, correction 
R3 §p> remaining loss coefficient 
R2 ; Oss coefficient due to wa 
We of Yriction 
RPRO Sp» total profile loss coefficient 
CL Rotor. tip..clearance 
¥CL Tip clearance loss coefficient 
RTOT Total-loss..coefficirent 
s Blade spacing 
DEZ Normal trailing edge thickness 
HM Blade height 
ak 3 ae 
PSID Se. esses due to fan 
PSIF ee 
UM Tip speed 
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TABLE A-X 


FORTRAN SYMBOLS IN SUBROUTINE FLOWR 








FORTRAN SYMBOL DESCRIPTION 

PRATGR Critical pressure ratio 

PHICR PORIT? Eritical Llow-function 
HSTAR H***, energy parameter 

XI Z, area reduction coefficient 
PHI 6, flow function (unchoked flow) 
ARAT Streamline throat DIA/mean throat DIA 
AS Mean stator throat diameter 

AR Mean rotor throat diameter 

WREQ M required to satisfy continuity 
WSUM M calculated 

WPER of M at each streamline 
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TABLE A-XI 


FORTRAN SYMBOLS IN SUBROUTINE ROTOL1 

















FORTRAN SYMBOL DESCRIPTION 
OMEG. w, wheel speed (RAD/sec) 
og oR tator mean 
U2 ut wroter mean 
“stator mean 
WU1 Wuyi see figure A-3 
BETA1 Br, see figure A-3 
W1 Wi, see figure A-3 
TTE Equivalent temperature 
PTE Equivalent pressure 
HE Equivalent enthalpy 
ZETAR Sp rotor loss™~coetfrcrent 
ZETAPR f, rotor profile Toss Coctricient 
DHEDX Enthalpy gradient between streamlines 
DSDX Entropy gradient between streamlines 
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FORTRAN SYMBOL 


TABLE A-XII 


FORTRAN SYMBOLS IN SUBROUTINE ROTO2 





DESCRIPTION 








BETA2 


8,, see figure A-3 


DBETDX dB 


Eee between adjacent streamlines 
Vw steXlal Velocity 

a2 
W2, see figure A-3 
Axial distance between stations 
Radial component of velocity 
Wus , see figure A-3 
Vu2, see figure A-3 
Relative Mach No.at rotor exit 
T2 
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2S 


P2 


[Total-static pressure ratio] 
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INPUT DATA AND APPROXIMATIONS 


CALCULATE OUTLET ANGLES 
COMPUTE M 


SOLVE EQUATION OF MOTION AND COMPUTE 
STATOR EXIT PROPERTIES 


COMPUTE STATOR LOSS COEFFICIENTS 


ADJUST M. 








COMPUTE ROTOR INLET CONDITIONS 


SOLVE EQUATION OF MOTION AND 
COMPUTE ROTOR EXIT PROPERTIES 








ADJUST M, 





RE-CYCLE ALL OF THE ABOVE 
ACCOUNTING FOR STREAMLINE CURVATURE 


PRESS. 
RATIO 
OK? 


YES 
PRINT RESULTS 


END 


FIGURE A-1: PROGRAM FLOWCHART 
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z CAMBER LINE LENGTH 
= CHORD 

2 THROAT DIAMETER 

z CURVATURE RADIUS 

2 SPACING 

= MAXIMUM THICKNESS 


TE= TRAILING EDGE THICKNESS PROJECTED IN 
PERIPHERAL DIRECTION 
TN= TRAILING EDGE THICKNESS, NORMAL TO FLOW 
DIRECTION : 
FIGURE A-Z: 


BLADE NOMENCLATURE 
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FIGURE A-3: VELOCITY DIAGRAM NOMENCLATURE 
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MAIN PROGRAM 





esse 


INTERCONNECTION OF THE SUBROUTINES IN THE 
TRAUPEL METHOD 


FIGURE A-4: 
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‘FIGURE A=5%  ,ENETIAL PROFILE -LOSS,.COEFFLCTENT AND MACH NUMBER 
CORRECTION FROM TRAUPEL 
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FIGURE A-6: T.E.THICKNESS CORRECTION FACTOR, MIXING LOSS 
COEFFICIENT AND FAN LOSS COEFFICIENT FROM 
TRAUPEL 
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FIGURE A-7: "REMAINING" LOSS COEFFICIENT FROM TRAUPEL 
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CALCULATE CRITICAL PRESSURE 
RATIO AND FLOW FUNCTION 
DU FOR EACH S.L. 









XE= ] - ie ca 


o) = CRITICAL FLOW FUNCTION 


ADD TIP CLEARANCE 


YES YES LOWER VA(3) 
VA(3) IS OK 


RETURN YES 


RAISE VA(3) 


RETURN 


FIGURE A-8: SUBROUTINE FLOWR FLOWCHART 
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CALCULATE ROTOR OUTLET ANGLES 
FOR STREAMLINES 2 AND 4 


CALCULATE = FOR EACH STREAMLINE 


DO FOR EACH S.L. 


SOLVE THE EQUATION OF MOTION 
AT THE ROTOR EXIT (STATION 2) 








NO 1 
ADJUST Y 
HAS 
SOLUTION 
CONV. 
2 
YES 
CALCULATE ROTOR EXIT CONDITIONS 
RETURN 


FIGURE A-9: SUBROUTINE ROTO2 FLOWCHART 
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FIGURE A-10: TIP LEAKAGE LOSS COEFFICIENT PLOT FROM TRAUPEL 
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APPENDIX: 05 


DERIVATION OF EQUATIONS USED IN THE PROGRAM 
B-1. EQUATION OF MOTION FOR RELATIVE FLOW: 
The -equation of motion for. relative flow Ref. [5] is 
——> —> —> 
VHe= WX (VxXW+aw) + Tw S 
(Set) 
Using cylindrical coordinates, the terms of EQN (B-1) may be 


expressed as follows: 





Beto, che > de 
Vie = mene as or (2) 
Lo L2 tr 
Fr r 
VXW = eee a ae poe 
= o2 or 
r Ma We Wr 


ig | OW Wa}, ce | (rw) _ OW, 
a2 or a or a6 


Ur o Wa = 2(r Wu) 
Fr 20 as (B-3) 
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wx (7 XW) = 


LO 





Ls ie 
ao We Wr 
WWr Wal |! | dCrWe) oWel | | We 2(rWu) 
d2 or ig or oo} + | 30 a2 


1 f oWe  dlrWw) 1 (30rWw) QWr \i + 
oo a (2 ~ 2a yt (EM eae 





[ve (Bie = BB) (RE ea + 
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1 f 9Cr Wu) a [2M _ aWa 
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(B- 5) 
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(B-6) 
Combining equations (B-{) through (B-6) the terms in (B-2) 


can be written as: 








(B-8) 


ns es 
SS 
— 





2 Wr aAWe si aS 
VW/ — ee a ————— SES 


(B-9) 
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Since the flow has been assumed to be axisymmetric, all deri- 
vatives with respect to 9 are zero. Thus, Equations (B-7), 


(B-8) and (B-9) reduce to, respectively: 











- 3 a an oe or (B-10) 
oHe = We Qo Wr o Wa War 
o> a2 oi .vet r 


——————— (Be 11) 


— ne —— EEE —= 





or (B-12) 


Equation (B-10) may be written as 


OCW) 2 Wr dre) og Wop Wr 

> 7 Wa or Wa; (B=13) 
Substituting into equation (B-11), 

due Wa Wr 2lrWa) yyy, Wr 

o% Fr We or oe 





\W o We s WwW Wue Wr * ois 
aie We a az (B-14) 
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Multiplying equation (B-9) by Wr and Equation (B-14) by Wa 
results in 


dHe _~ WueWr _9(rWu) Wr 
“SS ie cael 





WaWr 2We . au WiW. + W, T 2S 





ar “or! (B-15) 

and 
dHe —~. WaWr alrWw) 2\Wr 

CE aa es — oo 
W.W, —2We _ 2 wWaeeWe +t We T -2Se 

or a2 

(B-16) 
Adding these two equations yields 
dHe oHe as. 

Wr a + Wa cs | = = = 2+ Wa 32 a 


Since the flow has been assumed to be adiabatic, the total 


relative enthalpy, H is constant along a-streamline. Thus, 


jeg 


ans dHR 
S Wa 
V He OQ= =Wa>>- + r or es 


and re-arranging, 


a He Sos Wr dHe 


92 Wa or 
(B-19) 


From equation (B-19), éq. (5-17) can*be written as 
Dae) _ ltr aS 
az Wee ar 
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(B-20) 


_- 


Bet sat 99 bomaaan ns od eer vols | ats *% 
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Substituting Eq.s (B-19) and (B-20) into equation (B-15) gives 





_ We dHe ~ _ We Wr oC, Ww) W ¢Wr 
Wa or | ee We. or 2 
2 Wo. Wee Wi, Ww as 
W — Qw iS wie 
a or a Wa r ar cay 


Multiplying Equation (B-21) by — 4 yields 
R 


oHe ae NA 0 (Fr War) Wa oWr 





We = rQwWWetT 2S 


vs (B=72) 


hrs Expression 1s identical to equation (B-21) and is the 
equation which must be solved. It must be put into a form 
which can be solved by the computer. Re-writing equation 


(B-22) given that 


2 
Wa 2We 1 awa) 
er A Or 


Vg} Wr, 2Wa 20rWa) |, 
or Wo. “BS 4 ee y or 


HW Naw —2-2HR aT = 2 = 


yields 





(B-23) 
The relative enthalpy can be written 
2 2 2: eS 
x4e c 240 4 
y J qe de (B-24) 
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The equivalent enthalpy, defined in ref. [1] is 





2 ra S 
le gy aah cee eee ag 
E 2G Z2qe J 


(Bi-25) 


Hence, the relative enthalpy can be written as 


ae 
He = He - lst 
aa 8 (B-26) 


Also, the turbine outlet static temperature can be written 


as 





te = 8 Pac 


@a2z7) 


Substituting Eq. (B-26) and Eq. (B-27) into Eq. (B-23) and 


applying Eq. (B-21) to the rotor exit, gives 








2 WW. 
90Ve") oy, Wra ng Wun 3 (War) | 
oh coe ow caem 
v4 
ro wias 32 S| Hebe Ve He 
2 Tt He ie ae: Cp 
Wa? Sa =O 
oC i= 
P cas (B-28) 
Given the relationships: 2 Ww, 
Tae AH a 
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Equation (B-28) can be written as 
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pe 
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a Fr ee (B-29) 
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and substituting ee 2uWh, into equation (B-29) 
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Multiplying Eq. 


(B-30) 


(B-30) by () results in the dimension- 
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tess form of Equation (B-29): 
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Introducing the non-dimensional quantities 
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Equation (B-31) is written as 


foe We LY as* 
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ay x aX ae Was m 
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he fourth term of Eq. (B-35) is 
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(B-36) 
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(B-36) by ( oe ) and observing that 
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2 l 
(| + TAN 6 = ); 





cos* A 
| as* 











mB dacy*) | _ 
1a 2 errr? a ax 


) 
2Trng 28 4 2 cins2ey Um SB coSB | QUmUa costs 
o X % Win ¥ Wes. ae 


2 2 a 
ee p a aia “sate |B ae) 
<a ie A 





105 


2 : 
 potzbyw od aap | 2-8) aoktaupe 
i, = Pa a — . a a 


ta \* 
4~v"% “ \ 
- my Si on 4 

= «4 


- - re 4 
fv 
im “oe a 
ee Ms 
y» hc ei 
+k: 4 mee cette 


i 


+ 
i 





A a , 
) 4 (ot-®) spl antylg 


ra L ie a* 


4% 
‘ 





va 


a 
ae 





wer - 
oe 
on 


To account for streamline curvature the following terms are 
introduced: 


i? 


2 Ae \* 
. alee @ 


Cos?) = 


(B-38) 
where A, the angle between the axial and radial components of 
velocity at a point, is approximated as the average value be- 
tween two stations. 


Also, 





kK Ree ie we 
es Wa Ok 
(B-39) 
where 6r is the streamline shift throught the rotor defined 


as 6¥ = TpotoR — TROTOR 


OUTLET INLET (B- 40) 


Ssubstitutine Eqs. (B-38) and (B-39) into (B-37) yields 


el Gistersas «|. SR L7+ (AR) | ds? 
“a =—cos 6 (2K Vm 28) - ( 2 < 
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To obtain a dimensionless equation, the term 
Ci = AGeT 


is introduced into Eq. (B-41) giving 


Pee he preven ceive) Gn de [LTH Be chs 
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Equation (B-42) is the form of equation of motion used in the 


computer program. 


B-2. EQUATION OF MOTION FOR ABSOLUTE FLOW 


The equation of motion for absolute flow 


(B-43) 
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Differs from the equation of motion for relative flow 


VHezwx(vxw +aw )+TVS 
_ (B-44) 
only by the term Wx2W which is the Coroiolis acceleration. 
To obtain the programmed form of the equation of motion for 
the stator, the previous derivation is followed, but with 


Y= 0, He becomes H, W becomes V, and 8 becomes a. 


B-3 THE AREA RESTRICTION FACTOR Z 

The condition at the outlet of a blade row with boundary 
layers on both sides of the flow channel is shown in Fig. B-l. 
The flow is considered to be turbulent within the boundary 
layer while, outside the layer, the velocity of the flow is 
the theoretical velocity. Assuming a power-law velocity pro- 


file, the velocity may be written, 
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(B-45) 


The mass flow rate exiting the blade row can be expressed as 
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where O,. and Ven represent the ideal conditions for an 


isentropic expansion through the blade row to the discharge 
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pressure Pd, which is assumed to be constant across the blade 
Spacing. The discharge angle of the flow leaving the blade 


row is closely approximated by the expression [Ref. 1] 


aI OW. 
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cosKt (8-47) 


Inserting Eq. (B- 47) into (B- 46) and reducing yields 
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Equation (B- 49) can be written : 
e [- Xe 
ae = te Z 
- -51 
| TH | Xe ( = (B-51) 
Substicuting Eq. (B- 51) into (B- 45) yields 
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Using the displacement thickness given by 


§ *2 6. (40) (7 Al 


(B-53) 
the mass flow ratecan be written as 
* 
mee a aes 
Mm = Crx VTH By 
(B-54) 


The loss coefficient, expressed in terms of average kinetic 


energy lost is 
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(B-55) 
where E is the actual kinetic energy of the flow, given by 
Vv we 
TH 
E = PraVrm (a- 2S) +a (aa a 
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(B-56) 
Substituting Eq. (B-51) into (B-56) gives 
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The energy thickness is written as 
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(B-58) 
The loss coefficient can therefore be written as 
Ce 
é ! = OL 
= $* 
25 a B- 
[- 2 (B-59) 


The area restriction factor Z, is the fraction of the flow 
area through which the uniform theoretical velocity would 


produce the actual flow rate, thus 
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(B-60) 
Defining the energy parameter (a form factor) as 
wKR 
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(B-61) 
using Equations (B-59) and (B-61), Eq. (B-60) becomes 
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where Ep is the profile loss coefficient. 
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B-4. THE ENERGY PARAMETER, ees 
In Equations (B-53) and (B-58) the denominator of the 


integrand is expanded using the binomial theorem, so that 
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The integral of Equation (B-58) is now written as 


| ah . Pn m Mm 
f as \ hi’ or a eens " 


& |- a ae 


3 Am | 
Xe A Pe bs, adh eo 


which, on integration becomes 
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which can be written as 
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In a similar manner, 
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Substituting Eq. (B-67) and Eq-7 (B-O8)¥ into Eq. (8-61). the 


equation for H*** used in the computer program is obtained: 
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FIGURE B-1: BOUNDARY LAYER EFFECTS AT THE EXIT 
OF A BLADE ROW 
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APPENDIX: .C 
PROGRAM SEGMENTATION ON THE HP-1000 


Segmentation allows large programs to be run on the 
HP-1000. The program is divided by the programmer into a 
main program and several segments, which are stored on the 
disc. Each segment and the main program are then compiled 
and loaded. When the program is executed, the main program 
and its segments are called into memory individually, and 
only as they are needed for execution. In this manner, a 
program can run in a partition which is smaller than that 
program's total size. 

When the main program has performed all executable state- 
ments, the first segment is called into memory by an EXEC 
call. The system then loads that segment from the disc into 
a memory block following the end of the main program. The 
process is illustrated in Figure C-l1. Note; the main program 
plus the largest segment may not together exceed 29 k. Once 
a segment is in memory it can call another segment. 

When executing, any segment can call any subroutine which 
is attached to the main program. It was this feature which 
allowed the present program to be run. All subroutines were 
placed within the main program. In fact, the main program 
consisted of nineteen subroutines and functions. A segment 


may not return to the main program. Communication of data 
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between the main program and the segments is accomplished 
through a common block. 

The four segments of the present program are ''MAIN", 
Ponont ws "PART Z™) and “PART.3". The manner in which control 
is passed from the main program to the first segment and from 
the first segment to the second is as follows: 


BLOCK DATA 


END 

PROGRAM THESS 

DIMENSION INAM (3) 

DATA INAM /2HSH, 2HOR, 2HT / 


CALL EXEC (8, INAM) 

END 

PROGRAM SHORT (5) 

DIMENSION INAM (3) 

DIMENSION NAME (3) 

DATA INAM /2HSH, 2HOR, 2HT / 
DATA NAME /2HPA, 2HRT, 2H2 / 


CALL EXEC (8, NAME) 
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END 
PROGRAM PART 2 (5) 

DIMENSION NAME (3) 

DIMENSION NAMR (3) 

DATA NAME /2HPA, 2HRT, 2H2 / 
DATA NAMR /2HPA, 2HRT, 2H3 / 


CALL EXEC (8, NAMR) 

END 

PROGRAM PART 3 (5) 

DIMENSION NAME (3) 

DIMENSION NAMR (3) 

DATA NAME /2HPA, 2HRT, 2H2 / 
DATA NAMR /2HPA, 2HRT, 2H3 / 


END 
The "(5)" after the program name indicates that it is a 
program segment. Note the manner in which the program name 


is put into a data statement using the Hollerith notation. 
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SEGMENT 1 
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FIGURE C-1: PROGRAM SEGMENTATION-ILLUSTRATION OF THE MAIN 
PROGRAM CALLING A SEGMENT INTO LOGICAL MEMORY 
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APPENDIX: 'D 


RUNNING THE COMPUTER PROGRAM 


If the reader is unfamiliar with the HP-1000 Computer 
System, references [11] and [12] should be consulted before 


attempting to run the program. 


D-1. DATA INPUT 

Using the editor, input the following data into segment 
SsHORT™.. 

1. Turbine operating conditions: referring to Table 
A-III, type in appropriate data in lines 66 through 69 and 
74 through 78. 

2. Special input data/program control parameters: 
referring to Tables A-IV and A-V, type in appropriate data 
in lines 83 through 98. 

3. Turbine geometry: referring to Tables A-1 and 
moni. type in data for stator and rotor in lines 103 through 


186. 


B-Z COMPILING THE PROGRAM 
1. To compile the main program type: 
:RU,FIN4 ,MAIN:: 25,-,- 
Zo  sompaie the first segment : 


:RU,FTN4,SHORT::25,-,- 


EL 
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Ss. Gampile the second segment: 
RU PING PART 2: 325,~5- 
4. Compile the final segment: 


oR EP ENaPART 3:3 255 ~ 5° 


LOADING THE PROGRAM 
Type 
RU, LOADR 
Tap return key 
Will display 
LOADR: 
Type 
OP,LB 
Will display 
LOADR: 
Type 
>RE,3MAIN::25 
Will display 
LOADR: 
Type 
TRE, COHORT? 225 
Will display 
LOADR: 
Type 
SRE} sPARI 27225 
Will display 
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LOADR: 

Type 
VE eR ous 25 

Will display 
LOADR: 

Type 
> END 

After the end statement, the loader will display that the 


program is ready for execution. 


D-4 RUNNING THE PROGRAM 

Type 

RUN, THESS 

The program will be executed and no further action by 
the operator is required. The computed pressure ratio of 
each iteration of the outer loop of the program is displayed 
on the screen as it is calculated. The operator therefore 
has some idea where in the iteration process the computer 


program is executing. 
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APPENDIX: —E 


DISCREPANCIES IN MACCHI'S PROGRAM 

i. Main program, lines 21 and 22; the value, of ICL has 
not yet been read. 

2. Main program, lines 163-166; the Traupel method of 
calculating gas outlet angles does not take the Mach number 
into consideration. However, in lines 163-164, the program 
is attempting to draw a parabola through points which repre- 
sent outlet angle as a function of Mach number. 

3. Main program line 281; the calling of subroutine 
SLINE is questionable. Parameters are transferred to that 
subroutine, but many of them have not yet been defined 
(HE, DHEDX, WPER2, DSDX1). These undefined variables will 
be set equal to zero by the IBM 360 and 370 computers. Thus, 
in line 10 of subroutine SLINE, the value of DWDX will be 
zero and in line 17, division by zero will occur and the 
execution of the program should cease. 

4. Subroutine ROTORI lines 22 and 26; the stator radii 
are used in the calculation whereas the rotor radii should 
be used. 

§. Subroutine ASOSI, line 107; the correct Fortran code 
is 


ZEAAPS(1) = .5 * ZETAS(T) 
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6. subroutine ALOSZ, line 121; the correct Fortran code 
is 
ZETAPR(I) = .5 *ZETAPR(1) 
72" subroutine ALOSZ, lines 1235-126; the stator radii are 
used in the calculation whereas the rotor radii should be used. 
8. Subroutine ANGAIN, line 14; the correct Fortran code 
is 
AO = ATAN(1. -XCL/H*CH*COS(ANG1) /COS (ANG2) * 
TAN (ANG2) + XCL/H*CL*COS(ANG1)/COS (ANG2) * 
TAN (ANG1) 
Note: Since reference [2] was published, Professor Macchi's 
program has been further developed by Professor Macchi under 
private sponsorship [Ref. 13]. The new code however, is not 


generally available. 
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APPENDIX F 


_ COMPUTER OUTPUT 
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